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I. Introduction 

F ORCED periodic flows arise in a broad range of aerodynamic applications such as rotorcraft, turboma- 
chinery, and flapping wing configurations. Standard practice involves solving the unsteady flow equations 
forward in time until the initial transient exits the domain and a statistically stationary flow is achieved. It 
is often required to simulate through several periods to remove the initial transient making unsteady design 
optimization prohibitively expensive for most realistic problems. An effort to reduce the computational cost 
of these calculations led to the development of the Harmonic Balance method [1,2] which capitalizes on 
the periodic nature of the solution. The approach exploits the fact that forced temporally periodic flow, 
while varying in the time domain, is invariant in the frequency domain. Expanding the temporal variation 
at each spatial node into a Fourier series transforms the unsteady governing equations into a steady set of 
equations in integer harmonics that can be tackled with the acceleration techniques afforded to steady-state 
flow solvers. Other similar approaches, such as the Nonlinear Frequency Domain [3,4,5], Reduced Frequency 
[6] and Time-Spectral [7, 8, 9] methods, were developed shortly thereafter. Additionally, adjoint-based op- 
timization techniques can be applied [10,11] as well as frequency-adaptive methods [12,13,14] to provide 
even more flexibility to the method. The Fourier temporal basis functions imply spectral convergence as the 
number of harmonic modes, and correspondingly number of time samples, N , is increased. Some elect to 
solve the equations in the frequency domain directly, while others choose to transform the equations back 
into the time domain to simplify the process of adding this capability to existing solvers, but each harnesses 
the underlying steady solution in the frequency domain. These temporal projection methods will herein be 
collectively referred to as Time- Spectral methods. 

Time-Spectral methods have demonstrated marked success in reducing the computational costs associ- 
ated with simulating periodic forced flows, but have yet to be fully applied to overset or Cartesian solvers 
for arbitrary motion with dynamic hole-cutting. Overset and Cartesian grid methodologies are versatile 
techniques capable of handling complex geometry configurations in practical engineering applications, and 
the combination of the Time-Spectral approach with this general capability potentially provides an enabling 
new design and analysis tool. In an arbitrary moving-body scenario for these approaches, a Lagrangian body 
moves through a fixed Eulerian mesh and mesh points in the Eulerian mesh interior to the solid body are re- 
moved ( cut or blanked ), leaving a hole in the Eulerian mesh. During the dynamic motion some gridpoints in 
the domain are blanked and do not have a complete set of time-samples preventing a direct implementation 
of the Time-Spectral method. Murman[6] demonstrated the Time-Spectral approach for a Cartesian solver 
with a rigid domain motion, wherein the hole cutting remains constant. Similarly, Custer et al. [15, 16] 
used the NASA overset OVERFLOW solver and limited the amount of relative motion to ensure static 
hole-cutting and interpolation. Recently, Mavriplis and Mundis[17] demonstrated a qualitative method for 
applying the Time-Spectral approach to an unstructured overset solver for arbitrary motion. 
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The goal of the current work is to develop a robust and general method for handling arbitrary motion 
with the Time-Spectral approach within an overset or Cartesian mesh method, while still approaching the 
spectral convergence rate of the original Time-Spectral approach. The viscous OVERFLOW solver will be 
augmented with the new Time-Spectral algorithm and the capability of the method for benchmark problems 
in rotorcraft and turbomachinery will be demonstrated. 

This abstract begins with a brief synopsis of the Time-Spectral approach for overset grids and provides 
details of the current approach to allow for arbitrary motion. Model problem results in one and two dimen- 
sions are included to demonstrate the viability of the method and the convergence properties. Section IV 
briefly outlines the implementation into the OVERFLOW solver, and the abstract closes with a description 
of the benchmark test cases which will be included in the final paper. 

II. Extension of Time-Spectral Method for Overset Solution Procedures 

The fundamental issue of providing Time- Spectral capabilities to an overset solver concerns the nodes 
that dynamically move in and out of the physical domain as a function of the relative motion between the 
surface geometry and the background Eulerian grid(s). Projection into the frequency domain has infinite 
support in time requiring special treatment for these nodes. Figure 1 demonstrates the issue for a figurative 
case of a one-dimensional piston moving in a sinusoidal path. 


x(t) = cos (cut) 



a b c 


(a) Piston Motion 


u(t) 



Figure 1: Figurative piston trajectory and solution histories. The piston “cuts” nodes as it oscillates over 
the fictional grid in Fig. la. Shaded regions in lb represent blanked intervals through which the solution is 
undefined for nodes of the corresponding color. 


Fictional solutions at nodes a, b and c in Fig. la demonstrate dynamic hole cutting by the relative motion 
of the piston. The piston never reaches node a and thus its solution is continuous for the entire time-history 
of the period. The piston blanks node b briefly during the middle of the period and node c twice. Thus, 
node b has one non-periodic time interval associated with it while node c has two - one each represented by 
the solid and hashed lines. The shaded regions serve to highlight the time over which each node lies outside 
the physical domain with an undefined solution. Thus, special treatment is required for nodes b and c while 
the standard Time- Spectral approach can be applied to node a. 

For a brief review of the standard Time-Spectral approach consider a general time-dependent scalar ODE: 

— + R(u) = 0 (1) 

The continuous solution is projected (u = $ _1 u) into the frequency domain, 

* 7 k = ( 2 ) 

from data at equispaced time samples, At = tj = jAt. The derivative operator is applied in the 
frequency domain, i.e. -^Uk = iukiik, and the result is projected back into the time domain leading to 

<LD<l> -1 u + R(u) = 0 (3) 
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where V = 4>D4> _1 is the transformed differentiation operator* 1 . Direct assembly of the transformed differ- 
ential operator, D, is derived in [7,8]. Iterating Eq. 3 in pseudotime, r, provides a steady-state solution 
procedure to an underlying unsteady, yet periodic, problem, 

-j-u + 7>u + R(u) = 0 (4) 

dr 

The majority of spatial nodes in an overset or Cartesian simulation will be treated with this standard 
procedure, but a highly accurate and robust treatment is still needed for the complement of nodes that 
undergo dynamic blanking. A natural approach is to fill the blanked nodes with solution data via spatial 
averaging and proceed with the standard Time-Spectral method (cf. [17]). While attractive for its simplicity, 
this approach is inconsistent as it propagates non-physical information provided by an alternative governing 
equation (spatial smoothing governed by Laplace’s equation) into the physical domain via the infinite support 
of the spectral basis. This is demonstrated by solving the linear advection equation 


du du 

di +a frc =0 


(5) 


given u(x, 0) = sin (cjx) and u(0,t) = — sin(cuat) for x E [a, b\. The analytic solution is u(x,t) = sin(uj{x — 
at)). The “geometry” is modeled by a gap in x through which the solution and analytic residual are 
linearly interpolated at every step. The interpolation error across the spatial gap (Fig. 2a) at one time- 
sample corrupts the solution throughout the period (Fig. 2b) via the infinite support of the spectral basis. 
Additionally, spatial interpolation inhibits the expected temporal convergence as demonstrated in Fig. 2c. 

Instead of spatial averaging, the current method partitions the time history into intervals of contigu- 
ously unblanked time samples and thus abandons the global temporal support for nodes exhibiting dynamic 
blanking (See Fig. lb). This is consistent with the physics of disjoint domains separated by an impermeable 
boundary. It is important to note that these temporal partitions must still be evaluated on a uniform lattice 
to avoid manipulating the spatial operator. Thus, we seek a robust and efficient method for evaluating the 
temporal derivatives for non-periodic functions on a uniform lattice. 

Several approaches are evaluated to achieve this goal. One option would be to compute derivatives with 
finite- differences or a localized differentiable basis such as wavelets. However, compact bases offer limited 
accuracy [18] and wavelet differentiation requires special treatment at non-periodic interval boundaries [19]. 
Temporal extrapolation techniques are another option that fill artificial data in the undefined portion of the 
given period using data from the physical portion of the time signal. Fourier continuation which extrapolates 
a non-periodic function into a periodic function on a larger domain [20,21,22], and compressed sampling 
that employs L\ minimization to solve an underdetermined system [23], are both too costly and lack the 
requisite robustness. The current approach projects the solution of the unblanked nodes within a single 
partition onto a global basis spanning the same interval. The constraint of evenly-spaced temporal collocation 
points eliminates traditional optimally- accurate orthogonal polynomials (e.g. Chebyshev) due to Runge’s 
phenomenon at the endpoints. Least-squares projections of orthogonal polynomials are more stable but 
expensive, and result in an overdetermined system whose projection does not interpolate the solution data. 
Rational polynomials are offered as a viable alternative for highly accurate approximations on evenly-spaced 
grids [24] and will be employed in the current work as they posses the properties of discrete orthogonality 
and spectral-like approximation order. An alternate basis may prove more successful in the future so the 
current approach will be kept general enough to use any projection operator, provided the corresponding 
analytic differentiation operator is available. In this work, solutions across each non-periodic interval are 
projected onto a basis of rational polynomials and differentiated accordingly. Fully periodic intervals are 
still projected in the Fourier basis, resulting in a hybrid approach that employs the optimal available basis. 

Floater and Hormann [25] introduced a set of barycentric rational polynomials with high rates of conver- 
gence and no poles. Additional efforts to explore rational polynomials and their utility as a pseudospectral 
basis for spectral collocation methods include but are not limited to [26,27,28,29]. Tee and Trefethen [30] 

Alternatively, the residual operator can also be projected into the frequency domain, and the resulting equations solved 
iteratively in the frequency domain [3,6]. This approach takes advantage of the Fast Fourier Transform which is an N log N 
operation as opposed to the N 2 matrix operation in Eq. 3. However, for the small N (0(10 — 100)) typically encountered in 
practical 3D Time-Spectral applications the difference is minimal and the time domain formulation proves more flexible 
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(b) Exact and computed solutions 


10 ° 


10" 1 


10 ° 10 1 10 2 10 3 

log(AT) 

(c) Convergence with N 

Figure 2: Spatial smoothing Time-Spectral approach 
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offer the following definitions for the rational interpolant, r(x), in barycentric form: 


r(x) 


T - J£h - fk 

X — Xk J ^ 

k=0 


N 


£ 

k=0 


Wk 

X-X k 


and its corresponding differentiation operator, D, such that = Dr where = r(xk)' 


(6) 


Djk — 


W k 

Wj(Xj-Xk) 

N 


if j ^ k 


~ E D 3i if 3 = k 

i=0,i^k 


The interpolation is reformulated as a weighted sum of nodal basis functions with the weights equal to the 
function value at the nodes fk = f(xj . c ): 


N w k 

f(x) = 'Djfk(t>k{x), with 4> k {x) = (7) 

k = o w k 

k = 0 X ~ Xk 

Since barycentric rational polynomials interpolate the solution data, the transformation operator is the 
identity matrix, §jk = <fik{%j) = $jk, implying discrete orthogonality. It follows that the transformed 
differentiation operator is identical to D: V = = D. The weights, Wk , control the approximation 

and stability characteristics of the rational interpolant and are provided by Floater and Hormann [25] to 
avoid poles and provide an approximation order of d. 

Given N + 1 samples choose the desired approximation order d and define: 


It := {0, 1, ..., N — d} and J a := {i G It ’• ol — d < i < a} 


(8) 


Weights guaranteeing an absence of poles are defined by: 


u>k 


i-\-d 

(-!)*-“ e n 

i(zJ k 


1 

\x k - Xj\ 


(9) 


Readers are directed to the derivation in [25] for additional details. An example of three basis functions 
corresponding to the first, second and fourth nodes with approximation order, d = 3, over 8 sample points 
(N = 7) is provided in Fig. 3 with the values at the nodes highlighted with circles. Note the interpolation 
property of the basis, <pk {%j ) = djk- 



Figure 3: Rational polynomial basis functions 0o(^) 5 0i(a?) and 03 (x) corresponding to nodes 0, 1 and 3, 
respectively, for approximation order, d = 3, and 8 sampling points. 
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Approximation and differentiation performance of the rationals are compared to third-order finite differ- 
ences, the Fourier basis, least-squares approximation using orthogonal polynomials on evenly-spaced nodes 
(M < 2 y/~N) and Chebyshev polynomials on Chebyshev nodes in Fig. 4. An even-odd harmonic function 
and Runge’s function (designed to expose numerical instabilities) are differentiated by each of the aforemen- 
tioned methods, along with convergence of differentiation error, e = ||^ — 77>/ 1| 2 , with increasing N. Fourier 
differentiation is not applied to the aperiodic Runge’s function. The approximation order for the rationals is 
selected as d = min(|iiVj + 1, 10) for this case, but optimal selection is problem dependent and an area of 
continuing research [29] . The rational polynomials approach spectral convergence for both the periodic and 
non-periodic problems. 



log (AO 



(a) Analytic function and derivative for 
f(x) = cos (cut) + sin(u;t) 


(b) Convergence of differentiation operators on periodic function 



log(iV) 



(c) Analytic function and derivative for 

/(*) = T+k^ 


(d) Convergence of differentiation operators on Runge’s function 


Figure 4: Convergence of differentiated periodic and non-periodic functions 


III. Linear Advection Test Case using the Hybrid Time-Spectral Approach 

A linear advection equation is used to test the hybrid Time-Spectral strategy. The barycentric rational 
polynomial temporal differentiation operator, D r , will be used on the spatial nodes blanked over some 
duration of the period and the Fourier differentiation operator, T>f, is used on nodes with a periodic interval. 
Information propagates from left to right with wave speed, a = 1, and thus analytic boundary conditions are 
supplied at x = 0 and along the right edge of the blanked region. The domain is discretized with N x = 101 
points in space and N m 21 time samples. Flow is initialized uniformly to zero outside the boundary nodes 
and third-order finite differences are used to approximate spatial derivatives with results plotted in Fig. 5. 
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T) f T) r T) f 

< ► u{x,t) 



0 0.2 0.4 0.6 0.8 


X 

(a) Complete space-time solution 



X 

(b) Solution and blanked interval at fixed time t = 0.4286 

Figure 5: Solution for linear advection test problem using the hybrid Fourier-rational polynomial Time- 
Spectral approach. The rational polynomial differentiation operator, D r , is applied to nodes in the central 
region of the domain that undergo blanking and the Fourier differentiation operator, Vf, is applied elsewhere. 

The maximum error in the blanked simulation is actually lower than the unblanked case (7.717 x 10 -3 
vs 7.757 x 10 -3 ) due to the analytic solution boundary condition on the right side of the blanked region. 
The largest solution difference upwind of the analytic boundary was 7.861 x 10 -4 demonstrating that the 
Fourier-rational hybrid approach essentially produces the same result in the unblanked region for this simple 
test problem. The same approach will be used in the hybrid Time-Spectral OVERFLOW overset solver, as 
discussed in the following section. 

IV. Implementation of Time-Spectral Approach in OVERFLOW 

One of the strengths of the Time-Spectral approach is that the spatial operator remains unchanged, 
making the scheme straightforward to add and maintain in existing code. The current work leverages the 
OVERFLOW overset flow solver [31], which provides a mature code base for solving the three-dimensional 
compressible Reynolds-averaged Navier-Stokes (RANS) equations using a structured finite- difference scheme. 
OVERFLOW uses an approximate-factorization algorithm for the implicit solution of the spatial derivatives, 
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and the Time-Spectral operator is easily integrated with this approach by the addition of a similar factored 
operator for the time derivative. This avoids forming and solving a x Nq x N system of equations 
coupled across time and space. Writing the full Time-Spectral approximate-factorization difference equation 
in “delta” form gives, 

[I + AtA x ] [I + At A y ] [I + AtA z \ [I + ArA t ] AQ = -At (i?(Q n ) + 5(Q n )) (10) 

where R(Q n ) is the standard spatial residual operator, and S(Q n ) is the Time-Spectral differential operator 
applied to the state vector at time level n. A x , A y and A z are the linearized spatial operators corresponding 
to the x—, y— and z— directions and A t = V is the temporal differentiation operator. The Time-Spectral 
approach as applied to OVERFLOW retains the identical spatial residual operator and avoids modifying the 
existing implicit operator because they are applied one time-sample at a time. The required modifications 
are limited to an additional linear solve (of dimension N ) for the Time- Spectral factored operator and 
computation of the source term. In other words, time is treated as an additional spatial dimension when 
solving for the steady state in the combined space-time domain. 

The computational memory requirements for a Time-Spectral solution are greater than a standard un- 
steady case because the solution and source term must be stored for each time sample. However, careful 
implementation limits the storage requirements and allows for a significant number of temporal modes even 
in complex three-dimensional problems. Ignoring the trivial storage and processing associated with the hole 
cutting, the current implementation in OVERFLOW stores a copy of the state vector Q n and Time-Spectral 
source term S(Q n ) at each time level, along with an additional working copy of the grid and metrics. Thus, 
the total additional storage is 2N^NqN + 12 Ah), where N x and N represent the number of spatial and 
temporal modes, respectively, and Nq is the number of conserved solution variables. An additional cost 
of the approximate-factored matrix for the temporal derivative is also accounted for depending upon the 
relative size of N x and N. Using a target problem size of N% = 125M grid points, which is representative 
of a three-dimensional rotorcraft or turbomachinery simulation, and 125 8-core Intel Nehalem nodes of the 
NASA Pleiades supercomputer, Fig. 6 shows the variation in total memory with number of temporal modes 
for our target problem. As seen, even limiting to 125 compute nodes there is still sufficient memory for up 
to 200 temporal modes, and larger problems similarly benefit from the favorable scaling. 



Figure 6: Total Memory Storage (GB) versus N for 125 Intel Nehalem Nodes on the Pleiades supercomputer 

V. Final Paper 

As mentioned in the Introduction, the Time-Spectral approach within OVERFLOW has already been 
demonstrated by Custer et al. [15,16], and as such the current work to-date has focused on developing 
the necessary numerics of the Time-Spectral approach for disjoint domains using the rational polynomial 
bases. The use of the rational polynomial differential operator in an implicit two-factor scheme leads to 
a conditionally stable algorithm. When combined with the spatial operators in a multi-dimensional finite- 
difference scheme, as above, the numerical dissipation of the spatial operators can move all the unstable modes 
left of the imaginary axis. This is similar to the behavior when using the implicit Euler time integration in 
a two-factor scheme [32]. In practice the conditional stability of the rational polynomial operator does not 
appear more restrictive than the existing stability constraints of the spatial operator, however, a complete 
stability analysis of the rational polynomial scheme will be included in the full paper. 
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Three benchmark two-dimensional problems are targeted for inclusion in the full paper. First, unsteady 
and Time Spectral simulations on single and overset grids are computed for an inviscid transonic oscillat- 
ing NACA 0012 airfoil and compared to experimental data found in [33] followed by a periodic tri-rotor 
configuration (See Fig. 7) and a stator-rotor-stator configuration (Fig. 8). 



Figure 7: Two-dimensional tandem airfoil configuration. Flowfield is periodic in freestream direction to 
simulate a rotorcraft system in two- dimensions. 



Figure 8: Example stator-rotor- stator configuration where rotor moves down relative to the stators. 


These problems are representative of more complex three-dimensional applications, while still computa- 
tionally tractable for a rigorous analysis. Further, the rotor-stator configuration contains benchmark experi- 
mental data [34]. In particular we are interested in the ability of the Time-Spectral approach to capture the 
complex wake interactions, and associated broad-spectrum frequency response, when using hybrid-RANS or 
LES turbulence modeling approaches for these problems. These computational studies will be included in 
the full paper. 
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